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Abstract. We consider the geometric numerical integration of Hamiltonian systems subject to both equality and "hard" inequality 

constraints. As in the standard geometric integration setting, we target long-term structure preservation. We additionally, however, also 

consider invariant preservation over persistent, simultaneous and/or frequent boundary interactions. Appropriately formulating geometric 

methods to include such conditions has long-remained challenging due the inherent nonsmoothness they impose. To resolve these issues 

we thus focus both on symplectic-momentum preserving behavior and the preservation of additional structures, unique to the inequality 

constrained setting. Leveraging discrete variational techniques, we construct a family of geometric numerical integration methods that not 

only obtain the usual desirable properties of momentum preservation, approximate energy conservation and equality constraint preserva- 

^— H tion, but also enforce multiple simultaneous inequality constraints, obtain smooth unilateral motion along constraint boundaries and allow 

^— H for both nonsmooth and smooth boundary approach and exit trajectories. Numerical experiments are presented to illustrate the behavior of 

f^^ these methods on difficult test examples where both smooth and nonsmooth active constraint modes persist with high frequency. 

^ 1. Introduction. Recent attention has focused on geometric numerical integration methods that closely 

^ preserve invariants of continuous systems | Sanz-Serna and Calvo 1994[ Hairer et al.|2002 Leimkuhler and Re- 
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ich|2004|. In particular, symplectic-momentum preserving integration schemes have been successfully applied 



over a wide range of applications. In the unconstrained setting, these methods exactly preserve the symplectic 
form and momenta, and maintain good long term energetic behavior by approximately conserving energy (up 
to an additive constant) over long spans of simulation. 

While extensions of symplectic methods to equality constrained systems have been extensively investi- 
gated [Hairer et al.|2002 Leimkuhler and Reich|2004|, appropriately formulating symplectic methods to in- 
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r^ elude inequality constraints continues to remain a challenging problem that has received hmited attention in the 

"trt literature. The difficulty in treating these latter cases stems from the inherent nonsmoothness imposed by such 

conditions. Stewart] f2000|, in particular, notes that, despite good conservative properties in the smooth setting. 



symplectic methods do not necessarily translate easily to nonsmooth problems. Direct extensions of existing 
symplectic methods (e.g., implicit midpoint method, Newmark methods, etc.) loose their good energetic prop- 
^N erties and often generate surprisingly nondeterministic errors when subject to inequality constraints ||Stewart 

[20001 . 



-Y~, In this work we consider Hamiltonian systems subject to both equality and inequality constraints. As in the 

p^ standard geometric integration setting, we target long-term structure preservation. In the inequality constrained 

CN setting, however, we additionally consider structure preservation over persistent, simultaneous and/or frequent 

r- — . boundary interactions. 

f^ With these goals in mind we focus both on symplectic-momentum preserving behavior and the preserva- 

^^ tion of additional structures, unique to the inequality constrained setting. In particular, we develop geometric 

'^ numerical integration methods that not only obtain the usual desirable properties of momentum preservation, 

^ approximate energy conservation and equality constraint preservation, but also enforce many simultaneous in- 

equality constraints, obtain smooth unilateral motion along constraint boundaries and allow for both nonsmooth 
and smooth boundary approach and exit trajectories. 



1.1. Problem Statement. We consider the general case of a Hamiltonian system subject to generalized 
constraints that restrict system configurations, q, to an admissible set, A (i.e., q e A). In general, this admissible 
set can be both nonconvex and nonsmooth. 

More concretely, when the admissible set is given by a manifold, it will often be defined by a set of 
constraint equations of the form 

f(q) = (/o(q),...,/„(q))^ = 0. (1.1) 

For these cases the resulting holonomic Hamiltonian system reduces to a set of differential algebraic equations 
(DAE) for which a wide range of existing geometric approaches may be suitable jHairer et al.|2002[|Leimkuhler| 
|andReich|2004l . 

In the following, we consider generalizations where the admissible set. A, can always be described by 
both a set of constraint equalities (as above) and constraint inequalities, given by functions of the configuration 
variables, 

g(q) = (^o(q),...,^,„(q))^>0. (1.2) 



We start by considering the underlying, unconstrained system, with its natural Hamiltoniar|j and La- 
grangian functions. 
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We can then compactly include constraint forces by adding an additional potential term that, when extrem- 
ized, will penalize all nonadmissible configurations to ensure constraint compliance. In particular, to enforce 
so-called "hard" or "exact" constraints it is standard to consider the most extreme penalization available; we 
apply the extended value indicator function that infinitely penalizes noncompliance. For an arbitrary admissible 
set A, the indicator penalty function is given by 
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(1.4) 



The full, nonsmooth, inequality constrained, Hamiltonian and Lagrangian formulations for the system are then 
obtained by augmenting the unconstrained system's natural potential function, V, with the indicator function 



on the admissible set, I^ | Clarke 1983 1. 

The corresponding generalized equations of motion for the constrained system are then given by the Euler- 
Lagrange differential inclusion (DI), 



Mq + VV{q)e-dI^{q). 



(1.5) 



This inclusion, obtained from the apphcation of (9, the generalized gradient operatorrl to the nonsmooth penalty 
function [ Rockafellar and Wets|199 8~|, indicates that feasible constraint forces must lie in the inward pointing 
normal cone to the admissible set at q, given by — (9/^(q). On the interior of the feasible set the normal cone is 
trivially the origin. In this case the Euler-Lagrange DI simply reduces to the standard unconstrained equations 
of motion. Elsewhere, on the boundary, where constraint forces may be necessary, the normal cone is is given 
by a nontrivial, negative span of the active constraint gradients at q. See Figure [TT| 




Fig. 1.1. An example of an admissible set and a few evaluations of its associated normal cones. 



Our problem is then to derive geometric numerical schemes that integrate the generalized Euler-Lagrange 



system in (1.5i forward in time while preserving momenta and maintaining the good approximate energy 



preserving properties that are standard for symplectic integration schemes in the unconstrained setting | Hairer 



et al.|2002|. In the inequality constrained setting, however, there is additional critical behavior to consider: 



1. Smooth trajectories along the admissible set boundary. Consider, as a simple example, a mass- 



particle moving with a tangential velocity across a smooth region of the boundary (Figure 1.2 (a)). If 



a force is pushing the particle into the boundary, the particle's trajectory should remain smooth along 
the boundary. No upward bounce or other type of normal motion should result. Likewise, unless 
dissipative forces are additionally applied, constraints should not apply "boundary capturing"-type 
forces that slow the particle's tangential motion. Nevertheless, existing symplectic methods generally 



^Throughout, to simplify discussion and visualization, we consider the case of separable Hamiltonians with a flat metric on configu- 



ration 
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Note that here and in the following we reserve the d notation to exclusively designate the generalized gradient operator |Rockafellar| 
|andWets|1998) . 
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fail in one or both of these categories. This has important implications since, in the general setting, 
many important physical phenomena are covered by the abstraction of smooth boundary motion. Such 
smooth boundary systems include oscillatory elastica resting on foundations, twisted elastic knots 
(such as DNA filaments), mechanical and biological grasping models, masonry structures, granular 
flow, as well as many other boundary relative systems that are important across structural engineering, 
robotics, mechanics, biology and other critical appUcation areas. 
2. Nonsmooth exit and approach trajectories. Many trajectories in the hard constraint context can not 
be resolved by finite forces. For instance, consider a mass-particle approaching a smooth boundary 



in the normal direction (Figure 1.2 (b)). To resolve the impact, the mass-particle's trajectory must 
be reflected about the surface normal and thus we require a nonsmooth transition. More generally, 
whenever a trajectory is non-tangential to the admissible set boundary, an instantaneous nonsmooth 
jump in state is demanded. Nevertheless many symplectic-based approaches require, by construction, 
that constraints always be resolved by forces. Nonsmooth motion, however, is an intrinsic aspect 
of many important physical models applied to better understand the complex behaviors exhibited by 
mechanical colUsions, multi-impact systems, shock-wave propagation in granular media, molecular 
dynamics (using hard-sphere constraints), geophysical phenomena (e.g., earthquakes, iceberg-calving, 
etc.) as well as other difficult, nonlinear systems. As such, appropriately resolving nonsmooth motion 
is likewise critical. 
3. Smooth exit and approach trajectories. Consider once again the simple example of a mass-particle 



traveling smoothly along a boundary subject to a force pushing into the boundary (Figure 1.2 (c)). If 
the particle encounters a sufficiently large concavity, it must be able to "break contact" and leave the 
surface. Likewise, if a similar mass-particle is on a ballistic trajectory, that approaches the boundary 



tangentially (Figure 1.2 (d)), it must not bounce back off, but instead maintain a smooth on-surface 
trajectory. More generally, smooth trajectories, on the admissible set boundary, should not, in any way, 
be constrained to remain on the boundary, while initially off -boundary trajectories, that encounter the 
boundary tangentially, should remain smooth. 
Finally, we note that many fundamental applica- m»-_ * yl 

tions of inequality constrained systems (e.g., granu- \\\N.v«j>>;k^^^ \ / , ^ 

lar systems, biomechanical locomotion, structural en- ^ ^^ 

gineering) intrinsically combine both smooth and non- ^^v^..^,^. > ,,^ ^™„„™. ^ 

smooth modes and thus require the appropriate resolu- """^ 

tion of both modes within the same framework. P'°- ^ ■'^- Examples of (a) smooth-boundary, (h) nonsmooth. 

(c) smooth-exiting and (d) smooth-approaching trajectories. 



1.2. Overview. In the following sections we start by discussing related approaches (i 1.3 1 and then begin 
to incrementally construct our numerical integration methods, step by step. We first consider the discrete vari- 
ational perspective (Q and note that a nonsmooth, discrete Hamilton's Principle for unilateral constraints has 
not been previously considered. To address this gap we propose a Discrete Nonsmooth Hamilton's Principal 



that leads to a Discrete Euler-Lagrange Inclusion (DELI) formulation (S2.2i. We then observe that standard 
discrete variational structure is not sufficient, on its own, to define a well-posed integrator in this setting ({ 2.3 1. 
The proposed DELI formulation thus provides a framework for numerical integration in addition to which 
further structure is required to compose a fully specified integration method. To instantiate the first such DELI- 



based method we next consider the structure-preserving behavior of both smooth (j |3.4[ i and nonsmooth ({ 3.5 1 
time-continuous constrained trajectories. With these observations in place, we then return to the time-discrete 
problem (Q, and show that introducing time-discrete analogues of smooth (fpl and nonsmooth (f|6| bound- 
ary structures to DELI allows us to construct a corresponding, symplectic-momentum preserving, smooth- 
discrete integrator and a momentum preserving, nonsmooth-discrete integrator with observed good long-term 
energy behavior Finally, in f|7] we derive a discrete, generalized variational integrator (GVI) for generalized 
inequality-equality constrained systems that resolves both smooth and nonsmooth modes. Further numerical 
examples illustrating these methods are presented and discussed in ^ 

1.3. Related Work. To date, approaches for the structure preserving integration of such inequality con- 
strained, nonsmooth Hamiltonian systems can effectively be divided into two strategies: 

1.3.1. Direct-Substitution Methods. Energy-momentum and symplectic-momentum methods have been 
extended to the inequality case by the direct substitution of a nonsmooth constraint force [Laursen and Chawla] 
[T9971 IKane et al.|[T999l |Stewart|[2000l [Laursen and Lovel|2002l [Pandolfi et al1[2002l [Deuflhard et aU2(m\ 
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Fig. 1.3. Point-Mass Example: In these figures we examine the behavior of direct substitution methods on the simple example of 
a one-dimensional point-mass subject to gravity and a single ground plane constraint at zero. We plot position along the y-axis and time 
along the x-axis. All examples shown here begin with implicit midpoint as the base, unconstrained numerical method. In (a), starting with 
q{0) = 1, p(0) = 0, we plot the trajectories of the midpoint-constrained direct method in blue, the endpoint-constrained direct method in 
red, and the exact solution 's trajectory is plotted in green. Note, as well, that all trajectories are indentical until the first impact. In (b) 
we plot the the trajectories of the end-point constrained direct method for three small perturbations in initial position. Note that while the 
endpoint-constrained direct method generates fully dissipative solutions for all cases, the degree of dissipation varies in an energetically 
inconsistent fashion. 



Khenous et al. 2008 Krause and Walloth '2009! Doyen et al. 20111 . Following the time continuous case 



discussed above, the extended value Hamiltonian for inequality-constrained systems is naturally composed by 
concatenating the unconstrained system's potential with the extended value indicator function. Nonsmooth 
extensions of standard geometric methods are then composed by direct substitution of the obtained nonsmooth 
constraint force, — (9/^(q), as a standard force, handled as dictated by each numerical method. (See ^B.l for 
examples.) 

Unfortunately, as foreshadowed above, direct substitutions generally destroy energy conservation and gen- 
erate time-step dependent, nondeterministic behaviors when applied to symplectic methods. In particular, as 
we prove in Appendix B, all direct-substitution, nonsmooth, one-step, symplectic methods can not, in general, 
preserve, either approximately, nor exactly, the energy of the Hamiltonian system. In practice, this results in 
undesirable time-step-dependent and position-based restitution errors that can produce both energy growth and 
dissipation behaviors, while momentum similarly drifts. 

Energy-momentum based direct methods and a variety of extensions and stabilizations are similarly an 
active and promising area of research | Doyen et al.'201 Ij. Although comparable behaviors to direct symplectic 
schemes are currently observed in practice fDoyen et al.|2011 1, drift and stability properties of these methods 
remain under investigation. 

Example: A One Degree of Freedom Particle System. To illustrate some of these issues consider the 
following simple example of a one dimensional particle with mass, m~\, configuration, q e M, subject to a 
simple unilateral constraint g(q) = q > 0, under gravity such that y(q) — 9.8q. (Also see Figure 1.3 ) Stewart 



1 2000| uses this simple example to illustrate some of the potential issues involved in applying direct methods 
to existing symplectic integration methods. Poor energy conservation behaviors are observed for this case by 
the author for both the (symplectic) implicit midpoint rule [Hairer et al.|2002) and the IMEX Newmark scheme 
from Kane et al. 1 1999). In particular, 'Stewart 1 2000) notes that the energy of these systems is not conserved 
and that the effective coefficient of restitution varies with state for both methods. 

We start by dropping the particle from a height of q{Q) = 1. In Figure 1.3 a), we plot obtained position for 



a range of algorithms along the y-axis and time along the .^c-axis. In blue we plot the midpoint-constrained direct 

method; e.g., using constraint forces of the form —dlp^{- — j"^)- This generates, within the same simulation, 
effective coefficients of restitution that vary widely between dissipation and unstable energy-growth. In general 
these effects vary the with time-step size and relative position with respect to the constraint boundary. In green 
we plot the exact solution. Note that both standard collision integration methods (discussed below) and our 
new method proposed in SectionlTlgenerate trajectories that closely match this solution. 

Finally, we also plot the endpoint-constrained direct method, e.g., using constraint forces of the form 
— 5/y^(q'+' ), in red. Note that while the end-point method likewise does not preserve energy, it generates purely 
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dissipative energy errors. For this reason, end-point constraints and/or additional dissipative stabilizations 
have often been favored for direct methods in the literature |Kane et al.|1999j Pandolfi et al.[2002j ,Deuflhard| 
et al.||2008 Krause and Walloth 2009 Doyen et al. 2011| . The degree of dissipation, however, varies in a 



non-physical manner In particular, as with the mid-point constraint, restitution for the end-point constrained 
direct method is likewise dependent on time-step size and the relative position of the particle with respect to 
the constraint boundary. This is illustrated in Figure [T73| b), were we plot the the trajectory of the end-point 
constrained direct method over small changes in initial position. Note that while these three trajectories are all 
clearly dissipative, the degree of dissipation continues to vary in an energetically inconsistent fashion. 
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Fig. 1.4. Pogo-Stick Example: In these figures we plot the energy for a simple, two-dimensional, two-point, mass-spring system 
vertically oriented, subject to gravity, bouncing on a ground plane constraint. Both methods investigated here begin with Stormer-Verlet 
as the base, unconstrained numerical method. We set m= \ for both mass particles. The linear spring stiffness is set to k= 10, spring 
length is 1 = 5, and gravity is g = 9.8213. In (a) we plot the energy for the first 546 seconds of simulation time, stepping at h= 10^', for 
both the standard collision integration method in red and our generalized discrete variational integrator (GVI) in blue. After this step the 
energy in the collision integration blows up. In (b) we continue the GVI energy plot out for 100,000 seconds of simulation time. 



1.3.2. Collision Integrators. Alternately, collision integration methods have also been proposed | Stratt 



et al.'1981VMcNeil and Madden'1982t |Heyes|1982|[Suh et al.|1990||Houndonougbo et al.|2000||Fetecau et~ 
2003 ; Bond and Leimkuhler 2007 J in which unconstrained symplectic schemes interrupt fixed-size steps each 
time a unique constraint is identified as becoming active (i.e., about to be violated). For each such event, an 
impulsive constraint force (in the form of a reflection) is computed which conserves (depending on the scheme) 
either a continuous [Houndonougbo et al. 2000[ , discrete [Fetecau et al. 2003 1 or numerical Hamiltonian | Bond 
|and Leimkuhler|2007| . These methods thus aU effectively apply adaptive time-step subdividivision to resolve 
constraints. Given these variable-step modifications, a trajectory-shadowing, numerical Hamiltonian system 



can no longer be defined and thus energy conservation can not be expected in the long-term LHaireretal. 
[2002| . 

A novel stabilization scheme, based on the application of two complementary strategies, was recently pro- 
posed to improve the energetic behavior of collision integrators | Bond and Leimkuhler||200 7 1 . Stabilization 
is primarily obtained by applying reflections that preserve a second order approximation of the selected inte- 
gration scheme's (generally Stormer-Verlet's) numerical Hamiltonian. A secondary improvement is then also 
obtained by stepping to and away from constraint events using consistent, higher-order integration schemes. 
This latter modification, however, is motivated by applications in which constraint events are infrequent. The 
cost of higher-order integration is thus ameliorated by the assumption that the majority of time-steps are taken 
constraint free using the base, lower order method. In our development we now consider highly constrained 
systems where we expect most, if not all, steps to encounter inequality constraints. In this context, any potential 
stabilization advantage is then given almost exclusively by the higher-order energy reflection. We will examine 
this further in §6.1.1| 

While these developments can improve the performance of collision integrators, they still demand that 
only a single constraint be active at any time. Moreover, all collision integrators require expensive root-finding 
routines to determine each such constraint event. Many applications (e.g., surgical simulation, structural en- 
gineering, robotic grasping, granular flow, etc.), however, demand the simultaneous resolution of potentially 
large numbers of simultaneously active constraints. Individual treatment of constraints in these contexts can 
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lead to ordering bias, unacceptable stability issues, and generally have an especially large detrimental impact 
in implicit methods. Even when such ordering issues are acceptable, (e.g., small time-step, explicit methods) 
the individual time-ordered evaluation of all active constraint events will be, as |Cirak and West| ]2005[ note, 
computationally intractable for many simulation problems (consider, for example, Zeno's paradox). 

As a generalization it is reasonable, however, to expect that extensions of collision integration methods 
should be possible, even if not practical, which simultaneously resolve many active constraints. Indeed, one 
such possible generalization is suggested in'Fetecau et al. 12003]. Even in this paradigm, however, the inter- 



leaving of fixed-step symplectic -momentum methods with a variable-step, energy method still generates drift. 
To understand the practical implications of these issues in applying collision integration methods, consider 
the energy behaviors observed in Figure |1.4| ^a), above, as well as Figures |6.1[ a) and |8Tfb) and the related 
discussion in §6.1| and §8.0.1| 

2. Variational Integration and DELI Integrators. Despite the above negative results, nonsmoothly 
constrained Hamiltonian systems are invariant preserving in the usual senses [Kane et aL]|1999). Thus it is 



still desirable to find geometric integration methods that do not introduce numerical artifacts into simulations 
(dissipative or otherwise). One fruitful direction of research, that we do not pursue here, is the ongoing investi- 
gation in extending the discrete energy-momentum conserving framework to the nonsmooth setting. Instead we 
pursue the complementary research program of extending symplectic-momentum preserving schemes to resolve 
inequality constrained dynamics. In particular, observing the natural variational structure implied by inequality 
constraints [ Rockafellar and Wets|1998| , we revisit the application of Discrete Variational Integration | |Mars-| 



|den and West|200T | to the inequality constrained setting. 



To date, two approaches towards extending Variational Integrators to the nonsmooth, inequality con- 
strained setting have been considered: 

In the first approach, Nonsmooth Variational Integrators are generated fKane et aL]|1999[ Pandolfi et al. 



|2002[ by extremizing the standard discrete-action obtained directly from the quadrature of the underlying, 
smooth, unconstrained system. This obtains a standard Variational Integrator into which a nonsmooth con- 
straint force is then directly substituted. We observe that this approach then obtains nothing other than a 
direct-substitution symplectic method and is thus covered by Theorem B.l and subject to the same drift and 



stability issues as discussed above in ^1.3.1 Moreover, we observe that this approach violates the fundamental 



"discretize-first, extremize-second" mantra of discrete variational methodology. 



In the second approach, Fetecau et al. 1 2003 J observe that nonsmoothness in the variational formulation 



can be avoided altogether by adding an additional degree-of-freedom that parameterizes time of collision 
Composing and extremizing a corresponding discrete variational principle then generates Variational Collision 
Integrators which can be shown to be symplectic [Fetecau et al._2003J. As observed by [Bond and Leimkuhler 



1 2007), however, resulting schemes exactly obtain a standard collision integrator, as discussed above in n.3.2 



composed so that each collision preserves a discrete (rather than the continuous or numerical) energy at each 
reflection. 

With these observations in place we then note that, to the best of our knowledge, a nonsmooth, discrete 
Variational Principle for inequality constrained systems remains unexplored. We thus begin our development 
of nonsmooth geometric methods by developing these initial steps in the fully "discretize-first" manner. In 
particular, we apply Variational Integration (VI) methods [Marsden and West||2001[ [Hairer et"aLl|2002| as a 
prescriptive approach for generating a family of geometric methods for inequality constrained dynamics. 

2.1. VI Methods. Recall that VI methods, mimicking the continuous derivation of the Lagrangian equa- 
tions of motion, extremize a discrete action on the time interval t E [0, T], over all possible discrete paths of the 
form {q'^,...,q*^, ...jq'^}, where A^ = T /h. We begin with the discrete Lagrangian 

L,/(q^q^+')~^^ L(q,q)fl'f, (2.1) 

that gives a quadrature of the natural Lagrangian over a finite interval. 

The sum of the discrete Lagrangians, over all intervals, then gives the corresponding standard discrete 
action, 

'Y,LAq',q'+'). (2.2) 

*:=0 



Extremizing this associated discrete action then generates a unique, symplectic-momentum preserving integra- 



tion scheme | Hairer et al.[2002J, with accuracy specified by the order of the quadrature applied |West2004J. 



2.2. DELI. As a point of simple yet, to the best of our knowledge, novel departure we now add the 
nonsmooth penalty term to the discrete action and formulate a discrete, nonsmooth, constrained Hamilton's 
Principle that enforces hard inequality constraints on all endpoints of the discrete trajectory. 



N-l 



5,^ L,(q^q^+')-/A(q'•+') 9 0. (2.3) 



k=0 

Stationarity then gives us our constrained. Discrete Euler-Lagrange Inclusion (DELI), 

D2L,/(q'-',q')+£'li./(q',q'+')~5/A(q')3 0, (2.4) 

q'+' e A. (2.5) 

Given the last two configurations, q' and q'^', the DELI integrator advances state to a new configuration, q'+^ 
Here, in this recursive form, we observe that discrete constraint forces, applied to enforce constraints at the end 
of time-steps, are generated by normal cones constructed at the beginning of time-steps. 

We can further clarify this relationship by applying a momentum matching argument | West|2004) (or, if 



preferred, a discrete Legendre-Fenchel transform) to obtain the corresponding discrete momentum variables, 
p' and p'^^ as well as the discrete phase-space map composed of a position-map inclusion, 

DiL,/(q',q'+i) + p'-5/A(q')9 0, (2.6) 

e A, (2.7) 



^r+l 



and a corresponding momentum-map update equation, 

p'+'=D2L,(q',q'+'). (2.8) 

A one-step formulation of the DELI then follows. Given the discrete phase-space pair, (p',q'), at time 



f, we now apparently can solve the forward position map system, (2.6i and (2.7i, for the new constrained 
configuration q'+'. Then, with the new configuration in hand, the discrete momentum update follows directly 
from dOll. 



2.3. The Underdetermined Forward Map. We are now left with an interesting question of causality 
with respect to constraints. In our above DELI we require that corrective forces, that appear to be generated 
by the configuration at the beginning of the time-step (i.e, at q'), enforce constraints on the configuration at the 
end of the time step. 



In the direct methods we discussed above in Section 1.3.1 constraint forces are instead generated by 
configurations defined either somewhere in the middle of the time-step, i.e. q' + aq'+' , a £ (0, 1), or else at the 
end, i.e., q'+^. When solving inclusions with these terms, the variational structure of the generalized gradient 
then entirely specifies the constraint force. This, in turn, leads to the destruction of the invariant preserving 
properties of the underlying unconstrained method. See Appendix B, for further discussion. 

To retain the geometric properties of unconstrained methods we thus apparently require some additional 
degree(s) of freedom in choosing our constraint forces. Due to the odd time-causality noted above, the vari- 
ational structure of our DELI formulation no longer fully specifies constraint forces and thus provides such 
latitude for their selection. To see this, we first observe that in the DELI system merely requires constraint 
forces to lie in the span of the normal cone generated by the configuration at time t, i.e, — 5/y^(q'). Since this 



cone is independent of final configuration, (2.6 1 only specifies the span of possible constraint force directions 
but does not, unlike in direct-substitution methods, pin down a corresponding force magnitude. 

Thus, in the inequality constrained setting, we observe that standard discrete variational structure is not 
sufficient, on its own, to define a well-posed integrator. We find that constraint forces must be applied along 
the directions positively spanned by the normal cone and, likewise, all inequalities must be satisfied, but funda- 
mentally, nothing further is specified. In particular, an underdetermined system is composed with any number 
of solutions, most of which will not generate symplectic-momentum preserving maps. Our proposed DELI for- 
mulation thus provides a framework for numerical integration in addition to which further structure is required 
to compose a fully specified integration method. 

1 



In the following sections we will propose the first instantiation of a DELI-based method. We specifi- 
cally note that the dichotomy between nonsmooth and smooth trajectories is not well-resolved by the discrete 
variational framework. We will formulate our DELI-based integrator by first considering both the structure 
of smooth trajectories and the corresponding nonsmooth case. We will then show how imposing discrete 
analogues of these nonstandard structures upon the DELI formulation composes a fully defined DELI-based 
integration method and consider the behavior obtained by their application. 

3. Time-Continuous Setting. In order to consider how to design such methods, we first focus on the 
ways in which the flow of the constrained Hamiltonian system is transported in the time-continuous setting. 

3.1. Tangency and Normality. As a first step we introduce a few useful notational conventions and then 
briefly discuss variational structures that allow us to generalize notions of tangency and normality for cases 
involving multiple active constraints. 

3.1.1. Notation. For convenience, throughout the remaining sections we adopt a constraint subset nota- 
tion. Presuming, as above, that the a full constraint set is indexed by / e {0, ...,;«}, we then define, for each 
index subset, K C {0, ...,m}, the corresponding constraint subset, 

&^i^) = {\i^),-,S,,i^)y, {^o,..,^/}^K, (3.1) 

constraint subset gradient. 



G^(q) "=^' (Vg^^(q)^...,Vg^.^(q)^)^ {^o, ..,^/} ^ K, (3.2) 



and the active constraint subset gradient, 

N^(q) "^ G^j,(q), AK={/: g.{q)=0, ieK}. (3.3) 

Consistent with the notation introduced in Section [lT] for the vector- valued constraint function, g(q), an ab- 
sence of subscripting indicates that the entire constraint set is considered. We then have G(q) = G.^ , (q) and 

N(q) = N,„ ,„j(q). 

3.1.2. Normal and Tangent Cones. The nonsmooth generalization of a normal, the inward pointing 
normal cone | ,Rockafellar„1970; ,Clarke|1983[ [Rockafellar and Wets ,,1998 J to A at q, can, in our context, be 
defined directly with respect to the gradients of the constraint functions, g. It is given, with respect to the full 
active constraint gradient, N(q), as {N(q) x : x > 0}. Alternately, we note that for q in A, the generalized 
gradient of the indicator function generates the outward pointing normal cone to the feasible set. Thus, for all 
configurations we have 

^/A(q) = {-N(q)x:x>0}. (3.4) 

We note that, in particular, the normal cone reduces, on the interior of A, to the singleton, {0}. Correspondingly, 
recalling the Euler-Lagrange DI, given in \\.5\ , trajectories on the interior of the admissible set remain governed 
by the standard Euler-Lagrange equations of motion. 

Polar to the normal cone, the tangent cone to A at q, can also be defined with respect to the active constraint 
gradients, N, as 

T(q) ''J:' {y:N(q)^>0}. (3.5) 

When the admissible set is locally regular, the tangent cone, T(q), encodes the set of feasible directions at q, 
along which infinitesimal motion is locally permissible while, more generally, as we will discuss below, tangent 
cones generate useful generalizations of tangent spaces to nonsmooth regions of the boundary. 

3.2. Constrained Trajectories. In the inequality-constrained setting there is a critical dichotomy between 
nonsmooth and smooth motion. Consider the simple example of a mass-particle constrained to lie in a nons- 



mooth valley illustrated in Figure 3.1 (a). If the particle lies on the valley floor, then it could have arrived via 



some motion that is locally tangent to the valley's boundaries, such as moving along the valley floor, e.g., the 



green arrow in Figure 3.1 (a). In this case smooth trajectories, such as continuing along the valley floor, exist 
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Fig. 3.1. A simple motivating example of smooth and nonsmooth trajectories for constrained systems. See ^3.2\ and ^3.4\ f'or discussion. 

that can continue to satisfy admissibility. Alternately, the particle could have arrived at the valley floor with a 



downward pointing velocity, e.g., the blue arrows in Figure 3.1 (a) and in the cut-away view in Figure 3.1 (b). 
In this case there are no possible smooth trajectories that can continue to satisfy admissibility. Some sort of 
nonsmooth jump at the boundary is then required. 

To allow for the possibility of nonsmooth jumps we consider the left and right limits of velocity, 

q"(f)=limq(T) and q+(f) = limq(T), (2 6) 

Tp tit "- ^ 

where we assume that q has bounded variation on finite time intervals; limits on momentum follow similarly. 
We then consider trajectory intervals and events of interest with respect to the constraint set: 

Definition 1 . (Smooth Trajectory Intervals and Nonsmooth Trajectory Events). We identify smooth 
trajectory intervals as time intervals [a,b] C M where q(f) = q^(0 — (\^{t) far all t S [a,b] and nonsmooth 
trajectory events as times f G M where q^ (f ) ^ q^(0- 

3.3. Kinematic Conditions. Now consider the cut-away view of our point-mass example again, shown 



in Figure 3.1 (c) with both the tangent and negated tangent cones depicted in red. If the mass-particle lies on 
the boundary at time f , then it must have approached the boundary (in the left time limit) with a velocity in the 
negated tangent cone, since no other possible trajectory exists that could have reached the boundary. Likewise, 
in order to remain an admissible trajectory, the particle must leave the boundary (in the right time limit) with a 
velocity that is non-negative with respect to all active constraint gradients; thus the right-limit velocity should 
lie in the positive tangent cone. 

This intuition is made rigorous by a Lemma due to' Moreau]p988| Proposition 2.2], which says that a 
trajectory lies in the admissible region /or all time if and only if left velocities always lie in the negated tangent 
cone and right velocities lie in the positive tangent cone: 

Lemma 3.1. 

q(0 € A, far all t <(==> q"(f) G -T(q) and q+(f) € T(q), far all t. (3.7) 



For the nonsmooth case, as in the simple example discussed above, we expect that the left limit will not 
equal the right and thus we encounter a nonsmooth event requiring an instantaneous jump. We will discuss 



the time-continuous, nonsmooth case more thoroughly below in Section 3.5 First, however, we consider the 
time-continuous, smooth case in greater detail. 

3.4. Smootli Trajectories. We note that smooth motion, with respect to inequality constraints, corre- 
sponds at all times to either (unconstrained) motion on the interior of the admissible set, or to any combination 
of the following boundary cases: (1) encountering constraint boundaries tangentially, (2) leaving constraint 
boundaries, or (3) tangential motion along constraint boundaries. 

More concretely, smooth trajectory intervals are characterized by a tangent subspace condition: 

Lemma 3.2. If[a,b\ is a smooth trajectory interval then 

q-(f)=q+(f)e-T(q(f))nT(q(f)), far all t C,[a,b\. (3.8) 



Proof. Follows directly from Lemma 3.1 and Definition 1. D 

The intersection of the tangent cone and the negated tangent cone thus defines the tangent subspace at 
nonsmooth portions of the admissible set boundary (see our example Figure |3.1| (d)), along which smooth 
motion is possible and trivially generates the tangent space to smooth regions of the boundary. Conversely, we 
note that if the left velocity lies in the tangent subspace then, locally, an admissible smooth trajectory exists. 

In particular, while the stru cture of the normal cone implies that the standard (pointwise) Signorini- 
Fichera |Kikuchi and Qden|l988| complementaritjrl condition, g(q) _L A, holds everywhere, over smooth tra- 



jectory intervals Lemma 3.2 and the Euler-Lagrange DI in ( L5 1 imply that a stronger constraint-force continuity 
condition also holds per constraint: 

Theorem 3.3. Let gi[q{a)) ^QandVgi(^q{a)) q"(a) =0; ifXi{t) > holds for all t G [a,b] then 

(a) g,(q(f)) = 0, and 

(b) Vg,(q(f))^q+(f)=0, 
far all t € [a,b]. 

Proof. We can equivalently define the inward-pointing normal cone with respect to all constraint gradients 
as G(q)A — —dl^{q) with gi{q) > =^ A, — and A,- > =^ gi{<^) — 0. In turn, this implies that if A,(f ) > 
on an interval, i.e., t € [a,b], we must then correspondingly have gi{q{t)) = for all t G [a,b]. The chain rule 
then gives Vg,(q(f )) q+(q(f)) = on the same interval. D 

3.5. Nonsmooth Motion. With these observations on smooth motion in place, we turn to consider nons- 
mooth motion. As discussed above, when new constraints are encountered non-tangentially, i.e., from a normal 



direction, a nonsmooth trajectory event occurs (see Figure 3.1 (c)). In the general setting this corresponds to 
the case where the left limit velocity lies strictly in the interior of the negated tangent cone and requires the 
satisfaction of jump conditions: 

Theorem 3.4. Ifq'{t) ^ T(q) then 

(a) q+{t)^q-{t), 

(b) p+{t)^p-it)+dp, 

(c) dp e —dl^{q), and 

(d) q+(0€T(q). 

Proof. By Lemma 3.1 since q^(f) <^ T(q), q^(f) must lie strictly in the interior of — T(q) and, again by 
Lemma 3.1, right feasibility satisfaction of q+(f) e T(q) then requires a nonsmooth trajectory event at time 
t where q+(f) y^ q (0- To satisfy right feasibility, a jump, dp, in momentum is required to obtain p+(f) = 
p^(f) +dp £ MT(q). Finally, at t, the Euler-Lagrange DI reduces |Brogliato||l999) to the measure-valued. 



jump inclusion, dp e — (9/y^(q), which fully determines the feasible set of possible jumps. D 

At nonsmooth trajectory events, however, although the measure-valued Euler-Lagrange system, dp £ 
—dl^{q), continues to enforce momentum conservation, it no longer specifies energy behavior Conserva- 
tion of energy is then explicitly invoked by the side condition 

//(q,p+)=//(q,p-). (3.9) 

The combined criteria of the jump conditions and energy conservation are satisfied by all conservative, 
nonsmooth jumps at the boundary of the admissible set. In particular, these jump conditions uniquely define 
an energy and momentum preserving reflection whenever a single constraint is considered. In the general, 
multi-constraint setting, however, although the jump conditions effectively delimit the space of solutions, an 
invariant-preserving, constraint-enforcing, impulse is no longer uniquely given | |Glocker|2004) . In these cases, 
jumps can effectively be resolved by any of an infinite set of feasible corrective impulses. 

We note, however, that a geometric treatment of inequality constrained systems requires a well-posed 
resolution of jumps with respect to multiple simultaneously active constraints. Methods that satisfy well- 
posedness by providing uniquely defined, physically motivated, multi-impact solutions for the multi-constraint 
setting are thus desirable [Chatterjee and Ruina 1998J. In particular, in the following examples, we employ 



the generaUzed reflection operator of Kaufman et al. 1 2010 1. Alternately, the multi-impact solutions of Moreau 
1 1988) or Chatterjee and Ruina 1 1998| can also be applied. 



'Note that here and in the following we are using the standard complementarity notation, x _L y. to indicate xiyi = holds componen- 
twise for the matching entries of x, y G M"'^ ' . 
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4. Discrete Setting. Given the above analysis, we observe that nonsmoothly constrained Hamihonian 
systems correspond to smooth Hamiltonian systems almost everywhere. At instants of transition, where new 
constraints are encountered or left behind, either a smooth trajectory is maintained across the transition or 
a jump provides a nonsmooth mapping between two consistent smooth trajectories. For smooth transitions, 
motion is given by the Euler-Lagrange DI while, for nonsmooth transitions, jumps instantaneously reflect to a 



new, energetically consistent, symmetry preserving, smooth Hamiltonian flow. As discussed above in S 3.4 and 



\ 3.5 the conditions specifying the mode of each such transition are entirely given by the relationship between 
the left limit velocity and the tangent cone. 

To carry these observations into the discrete setting, we will thus first consider compatible, discrete- 
predicates for distinguishing between smooth and nonsmooth transition modes at each time-step. With these 
predicates in place, we will then generate compatible discrete analogues of the time-continuous, smooth and 
nonsmooth maps. 

As a prescriptive guide for constructing these maps, we next observe that, if we require our discrete maps to 
reduce to symplectic steps in the unconstrained case, any obtained unconstrained step ending on the constraint 
boundary is always trivially shadowed by an unconstrained numerical Hamiltonian system | Hairer et al.|2002). 



If the successive constrained time-step is discrete-smooth, then we should require the corresponding discrete- 
smooth map to be symplectic, so that it is likewise consistently shadowed. Similarly, if a successive constrained 
time-step is discrete-nonsmooth, then the corresponding discrete-nonsmooth jump-map should be constructed 
so that it transitions to a new, uniquely defined, invariant preserving numerical Hamiltonian system that is 
similarly consistent. Following this guide, we conjecture that a piecewise-smooth, numerical Hamiltonian 
system, for the full inequality constrained system, can then be constructed by gluing together the corresponding 
smooth numerical Hamiltonians over successive smooth and nonsmooth transitions. 

Below we will employ the preceding time-continuous analysis to show one way in which such discrete 
maps can be generated. The resulting integrator will, by construction, attempt to preserve correspondence to 
such a composite, shadowing Hamiltonian system. We note that one fundamental requirement of maintaining 
such a correspondence is that all new constraint boundaries be considered, in the discrete setting, only at the end 
(and/or beginning) of each time-step. Indeed, it is specifically the explicit violation of this particular criterion 
that causes collision integrators to perform so poorly in the nonsmooth casaj With these observations in place 
we now begin the development of our proposed, DELI-based, nonsmooth geometric methods. 

5. Discrete Smooth Setting. In the discrete setting, smooth-trajectory intervals are determined by consid- 
ering whether the discrete-analogue of the tangent subspace condition in ( |3.8[ ) is satisfied at interval endpoints 
along the discrete path such that 

M-ip'e-T(q')nT(q'). (5.1) 

We then identify each DELI forward-map with a discrete-smooth interval on the discrete path. Making all 
constraint functions and their gradients explicit in ( |5.1[ l we find that the discrete left-limit 

gM) = ^=^ Vgi{q'fU-'p' = 0, (5.2) 

and discrete right-limit 

g,(q'+')=O^^Vg,(q'+')^M-'p'+' =0, (5.3) 

subspace tangency conditions simply require momentum to lie in the cotangent spaces of all active constraints 
at the beginning and end of all smooth intervals. 

The discrete-analogue of the constraint-force continuity condition in Theorem 3.3 is then similarly given 
by 

A,>0=^g,(q'+')=0. (5.4) 



5.1. Discrete-Smooth Integrator. Reconsidering our derived DELI in ( 2.6 1 we first make constraint func- 
tions and their corresponding gradients explicit to obtain the equivalent forward map system, 

DiLrf(qSq'+i) + p' + N(q')A=0, A > 0, 

g(q'+')>0, (5.5) 

p'+i-D2L,(q',q'+l). 



''Of course, by construction, collision integrators are explicitly designed only for the nonsmooth case. 

11 




100 200 300 400 500 600 700 800 900 1000 



(a) 



(b) 



(c) 



Fig . 5 . 1 . As illustrated in (a), we consider a smooth trajectory of a constrained system composed of a simple, four-dimensional, two- 
point, mass-spring system, subject to a circle inequality constraint and an outward-pulling, radial potential. We initialize both particles 
to a rotational trajectory, initially tangent to the boundary and note that all constraints and potentials, in this example, are rotationally 
invariant and thus angular momentum should be preserved. Both methods investigated here begin with implicit-midpoint as the base, 
unconstrained numerical method and are stepped at h = 5x 10^'. In (b) we plot the energy (in green) and the angular momentum (in red) 
of the direct, end-point constrained method for the first 100 units of simulation time. We note the decay in both angular momentum and 
energy. In (c) we similarly plot the (approximate) preservation of energy (in green) and (exact) conservation of angular momentum (in 
red) obtained by our Discrete Smooth Integrator out to 1000 imits of simulation time. See ^5.2\below for further discussion. 



We then augment DELI with discrete left-limit tangency as a precondition, discrete right-limit tangency 
as a postcondition and discrete constraint-force continuity over each smooth interval to obtain a DELI-based, 
Discrete-Smooth Integrator. 



DiL^(q',q'+i) + p' + N(q')A-0, 
0<A±g(q'+i)>0, 



p'+i=D2L^(q',q'+'] 



N(q'+')M, 



N(q 



l+l\T 



M 



-Ipr+l 



0. 



(5.6) 
(5.7) 
(5.8) 
(5.9) 



An update step for t he Discrete-Smooth Integrator is then generated by first solving the positio n up date equa- 
tions, (5.6 1 and (5.7 1, for q'+' and X. A momentum update then follows from solving equations (5.8 i and (5.9 1 



(5.9) is symplectic-momentum conserving. 



for p'+ and /z. 

We then observe that 

Theorem 5.1. The Discrete Smooth Integrator given in ^5.6\ 
Proof. See fjC] 

By construction, the method likewise preserves the subspace tangency of discrete smooth trajectories when 
moving along, leaving and approaching constraint boundaries. Likewise, the integrator trivially reduces to the 
unconstrained symplectic method, given by the discrete Lagrangian, on the interior of the admissible set. 

5.2. Spring and Sphere Example. As an initial consideration of the Discrete Smooth Integrator and 



in comparison to the direct-substitution schemes discussed in ^ L3.1 we examine the behavior of a simple, 
smooth trajectory system. We consider the simple example of a four-degree of freedom, mass-spring system. 
Configuration is given by the positions of two mass particles, q ~ (qi,q2) G M^, connected by a spring. The 
system is subject to a standard spring potential between particles and an outward pulling, radial potential, i.e., 
y (q) = ^4=1 25/(qf q,) + j(|| qi — q2 || ~l)'^, with I — 2-^2. In addition, we apply inequality constraints that 
require the particles to lie in a sphere, g,(q,) =|| q,- || — r > 0, of radius r = 5. See Figure 5.1 (a). Note that 



all constraints and potentials, in this example, are rotationally invariant and thus angular momentum should be 
preserved by the system. 

We initialize the example system to a smooth, rotational trajectory: both particles start, at time f = 0, 
on the sphere boundary, qi (0) = (4, —3), q2(0) = (3, —4), with corresponding unit length momenta that point 
counter-clockwise and tangent to the boundary. This initiates an on-boundary, counter-clockwise rotation of the 
particles, as well as an oscillatory motion between them. Both methods investigated here begin with implicit- 
midpoint as the base, unconstrained numerical method and are stepped at /i = 5 x 10^'. 



In Figure 5.1 (b) we plot the energy (in green) and the angular momentum (in red) of the direct, end-point 
constrained method for the first 100 units of simulation time. We note that the decay in both angular momentum 
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and energy for the direct method corresponds to the destruction of both the rotational and oscillatory modes of 
the system. 



In Figure 5.1 (c) we similarly plot the (approximate) preservation of energy, in green, and the (exact) 
conservation of angular momentum, in red, obtained by our Discrete Smooth Integrator out to 1000 units of 
simulation time. Throughout this simulation both particles maintain the smooth, on boundary, rotational and 
oscillatory modes of the system. 

Finally we note that collision integrators are not suitable for application in such examples since they are, 
as discussed above, formulated to resolve only nonsmooth boundary interactions. Nevertheless if we do apply 
a collision integration scheme in this smooth setting (effectively applying a micro-collision model | |Brogliato| 



|19991), we observe that while energy behavior is not compelling (see E6.1.1 below for a similar example) the 



trajectory is also qualitatively wrong; the nonsmooth reflections applied by the collision integrator introduce an 
artificial, nonsmooth, normal, bouncing mode (off the sphere boundary) that increases over time and eventually 
drives the system to instability. 

6. Discrete Nonmooth Setting. If the discrete left smoothness predicate is not satisfied at time t, the left 
limit velocity must lie strictly in the negated tangent cone such that M^'p' ^ T(q'). A nonsmooth jump is 



then required to satisfy admissibility. Following our discussion in S 3.5 we apply discrete, generalized, energy 



preserving reflections to resolve these cases. Later, in ^6.3 we will focus on how we derive such discrete 
generalized reflection operators. First, however, we will consider the critical issue of when such reflections 
should be applied. 

6.1. Extended Reflections. We will, in this section, temporarily ignore smooth motion and presume that 
all encounters with the admissible set boundary are nonsmooth and are always resolved by the application 



of some energy preserving, generalized reflection (to be defined below in ^6.3i. Presuming that we apply 
an integrator that reduces to a symplectic method in the unconstrained case, standard backwards error anal- 
ysis [Hairer et al.||2002| guarantees that, away from boundaries, there always exists a trajectory shadowing 



numerical Hamiltonian. In an ideal setting each such unconstrained, symplectic time-step would end either in 
the interior or exactly on the boundary of the admissible set. 

Correspondence to the piecewise-smooth, numerical Hamiltonian, discussed above in fHj then follows. Of 
course, such nice behavior can not be expected in our generally unsynchronized universe. Instead, we are faced 
with the likelihood that most, if not all, new constraint boundaries will be encountered midstep. 



The standard approach of existing collision integration schemes (see S 1.3.2 1 is to take unconstrained sym- 
plectic steps until the admissible boundary is crossed. When a boundary is crossed, collision integrators return 
to the beginning of the time-step and then repeat the forward step using just the time-step fraction that lands 
configuration exactly on the active boundary. A energy-preserving reflection is then performed using current 
state and then, finally, the integrator completes the remaining portion of the full time-step. 

A simple alternate approach, that, to the best of our knowledge, we introduce for the first time here, 
is to always take full time-steps, using augmented or extended reflections to incorporate otherwise missed 
constraint events. More concretely, we apply smooth symplectic steps whenever the admissible boundary 
is not encountered from a normal approach. Anytime current configuration, q', indicates that a constraint 
is active and/or a predictor configuration, q'', indicates that a constraint violation is imminent or alternately, 
depending on the choice of predictor, has already occurred, we perform a generalized reflection at the beginning 
(equivalently end) of the time step, with respect to both active constraints and predicted active constraints. Once 
the reflection is applied we then take the full smooth, symplectic step. This ensures correspondence between 
numerical Hamiltonians (as discussed above in f|4]), at the cost of inducing errors in the representation of the 
local constraint boundary geometry. 

6.1.1. Pogo-Stick Example. To illustrate our new approach and to compare it with collision integration 
schemes, we consider the simple example of a two-degree of freedom, mass-spring system, subject to a single 



ground constraint. See Figure 6.1 (b). We apply Stormer-Verlet as the base, unconstrained numerical method. 
Configuration is given by the positions of two mass particles, vertically oriented, such that the positions of the 
top and bottom mass particles are given by q = {q\,q2Y G I^^- The system is then subject to gravity and a 
ground-plane constraint on the bottom particle, g(q) = ^2 > 0. We set m = 1 for both mass particles. The linear 
spring stiffness is set to fc = 10, spring length is Z = 5, gravity is g = 9.8 and the base, unconstrained step-size 
is/i=10-'. 



In Figure 1 .4 a) we plot in red the energy of the standard collision integration approach and in blue we 
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Fig . 6 . 1 . In (a) we plot the energy for the same simple, two-dimensional, two-point, mass-spring system, shown in (h), and discussed 
in Figure \L4\ and in ^6.1.1\ below. Here we use exactly the same physical parameters and algorithms as in the original example, with the 
exception that reflections for both the collision integrator and our extended reflection method are formulated to conserve Verlet's second- 
order numerical Hamiltonian, rather than the standard continuous Hamiltonian. We plot the energy for the first 1000 units of simulation 
time for both the collision integrator, in red, and our extended reflection scheme, in blue. Note that, after t = 713, the energy in the collision 
integration system blows up. 



plot the energy of our new extended reflection approach using a simple forward time-step predictor. In this first 
example both methods apply reflections that conserve the continuous Hamiltonian. The collision integration 
method, however, applies reflections at times of contact with the floor, while, in this simple example, our 
extended reflection method reduces to simply applying reflections at the beginning of all time-steps that would 
otherwise result in the bottom mass-particle passing through the ground-plane. 



In the first time-period, shown in Figure 1 .4 a), energy remains fully bounded for the extended reflection 
method proposed here and, in particular, varies in the expected manner of standard symplectic methods. Over 
the same period the hybrid solution experiences energy growth. After t = 546, energy in the collision integrator 



enters a highly oscillatory regime causing the solution to blow up. In Figure 1 .4 b) we continue the energy plot 
of our extended reflection approach out to 100,000 units of simulation time. 

Of course we are not limited to continuous Hamiltonian conserving reflections. We may, for instance, 
expect an improvement in behavior by applying reflections that more closely conserve the shadowing numerical 



Hamiltonian. In Figure 6.1 (a) we plot the energy of exactly the same two simulations discussed above except. 



in this example, we apply reflections that preserve Stormer- Verlet's second-order numerical Hamiltonian | Bond 



and Leimkuhler 2007] . We observe a small improvement in the bounds of the energy envelope for our extended 
reflection method. Similarly, we note that, although the collision integration scheme still experiences energy 
growth, the obtained solution does not blow up until later, at f = 713. Note that we could also consider applying 



higher order symplectic Gauss methods for unconstrained steps leading to collisions [Bond and Leimkuhler 
2007 1 . Given the persistence of constraint events, however, this is effectively the same as simply switching to 
a higher-order, base, unconstrained integrator for both methods. 



6.2. Discrete Limit Momenta. We next observe that applying the constraint force N(q')A in (5.5 1 with 
respect to the next time increment, f + 1, is indistinguishable from directly modifying the discrete momentum, 
p', at time t. We thus decompose the constraint multipliers X into right components, A+ — A' , given by a 
previous nonsmooth reflection, and left components, A^ = 1^'+'' , given by the current discrete-smooth step. 
Applying an energy conserving, generalized reflection operator at the end of time-step f , we then correspond- 
ingly obtain discrete analogues of the left and right limit momenta, 



P =P , 

p'^=p' + N(qOA + . 



(6.1) 
(6.2) 



Here N(q')A+ gives the impulse generated by the reflection of p' . Consistent with our discussion in ^6.1 the 
above partition of A requires generalized reflections to be applied exclusively at the beginning (or equivalently 
end) of each time-step. Given our above analysis, we then include the extended reflection in our definition of 
left and right discrete momenta. 

At the beginning of each time step, we augment the active set with the set of all constraints that would 
be activated or violated by a predicted configuration. More concretely, we predefine a consistent predictor 
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configuration, q'', and then define an extended active set. 

def 



(q',q'') 



{/: g,(q'')<0, /e{0,..,m}}u{j: g,(q')=0, ; G {0, ..,m}|. 
The corresponding extended discrete left and right momenta are then given by 



P =P 
P =P 



-GA(q')A^ 



(6.3) 

(6.4) 
(6.5) 



6.3. Discrete Generalized Reflections. To compute the discrete reflection term given in (6.5 1 above, we 
require a discrete generalized reflection operator, R, of the form 



p- + GK(q)A=R(q,p-,K,£), 



(6.6) 



that defines a unique reflection with respect to an active constraint set K. We require that any such reflec- 
tion operator satisfy the discrete analogues of the jump conditions from { 3.5 given by the discrete kinematic 
feasibility condition. 



GK(q)^M-'(p-+GK(q)A)>0, 
the discrete Euler-Lagrange inclusion condition, 

p+-p" = GK(q)A, A>0. 



(6.7) 



(6.8) 



and conservation of a specified energy function, E, (e.g., the discrete, numerical, or continuous Hamiltonian) 
such that 



Eiq,P^)=E{q,p ). 



(6.9) 



As in the time-continuous case, whenever multiple active constraints are considered, i.e., |K| > 1, this 
problem is underdetermined. From the infinitely many possible solutions, we select the one obtained by ap- 
plying a discrete adaptation of the generalized reflection operator of Kaufman et al. | 2010| . This operator 
simultaneously guarantees uniqueness, determinism, symmetry-preservation, and satisfies all jump conditions 
when multiple constraints are active. In SJAlwe briefly derive this discrete extension and provide pseudo-code 
for the implementation we employ. 

7. Generalized Discrete Setting. Returning to the full setting where both discrete-smooth and discrete- 
nonsmooth modes are possible, we finally consider a fully general DELI-based system. Starting with the 
DELI system given by (2.6 1, (2.7 1 and (2.8 i, recall that in f|5] we presumed, a priori, the precondition of the 
left discrete-smoothness criteria to obtain the Discrete-Smooth Integrator given by (5.6 1 through (5.9 1. In the 
general setting, however, at any time a subset of the active constraints may be discrete-smooth, while their 
complement may be discrete-nonsmooth. In such cases we can expect that momentum will be negative with 
respect to some constraint gradients, and thus nonsmooth modes will be present. 

To partition out constraint subsets that satisfy the discrete-smooth precondition at time t, we define the 
smooth constraint set as the set of discrete-smooth constraints, with respect to the discrete-right time-limit at t. 



§(q',p'^) '=1' {/: Vg,(q')^M-'p'"-Oandg,(q')-0, /e{0,..,m}}. 



(7.1) 



Then, adding the discrete left and right jump terms given by ( |6.4| l and ( |6.5| l, we include reflections and obtain 
a full, discrete, generalized extension of DELI for inequality-constrained systems, 



p' =p' + GA(q')A + , 



DiL,/(q',q'+') + p' +Ns(q')A-=0, 



0<A-±gs(q'+^)>0, 






p-'-=D24/(q',q'+') + Ns(q' 



;+h 



M, 



N§(q'+')'M-^p'+'=0. 
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(7.2) 
(7.3) 
(7.4) 
(7.5) 
(7.6) 



Algorithm 1 GVI(q, p,E,q'') II input: q', p', energy function, position predictor 

1: q"'^' -i— q // cache current position 

2: A<— j/: ^,(q'0<0, i ^{0,..,m}>\j\i : gj{q)=Q, _/' e {0, ..,m} ^ // compute extended active set 

3: p -^ R(q,p, A./i) // apply generalized reflection 

4:S-s-|/: Vg,(q)^M"'p = and §;(q) =0, / G {0,..,m}l. // compute smooth set 

5: q<— -^x: DiL(/(q,x) + p+Ns(q)A = 0, < A _L gs(x) > > // apply position update 

6:p<-|y: y =D2Lj(q"''',q) + Ns(q) ;U, Ns(q)^M"'y = o| // apply momentum update 

7: return q,p // output: q'+^ p'+' 



7.1. Algorithm. With the above developments in place, our proposed discrete, generalized variational 
integrator (GVI) for fully generalized inequality constrained systems is then given by the pseudocode below in 
Algorithm 1 . A GVI takes as input the current state, the selected energy function and the position predictor. It 
then computes the GVI step and outputs the next state. 

7.2. Bilateral Constraints. Finally, we now further presume that, in addition to the set of inequality 
constraints, we also wish to enforce a set of equality constraints, as given by \\.\\ . When we consider equality 
constraints, we observe that the admissible set is now given by a projection of the inequality constraints onto 
the manifold defined by the additional equality constraints. In turn, this implies that normal cones are now 
given by the negative span of the active inequality constraint gradients, projected down to the cotangent spaces 
of the equality-constraint manifold. 

It then remains to note that, in the generalized constraint setting, we treat equality constraints as active at 
all times. We then can add the equality constraints, f, and their corresponding gradients, which we will denote 
by F, to our discrete-smooth formulation, during all time-steps. 

For our GVI methods this requires three basic changes. First, for arbitrary inequality constraint subsets, 
K, the corresponding active inequality gradient subset must be projected onto the local manifold's cotangent 
space at q. Thus, when equality constraints are applied, we projecrlthe active constraint gradient and redefine 
as 



NK(q) - (/-(F(q)^M-iF(q)) + F(q)F(q)^M-i) GAK(q), AK={/: g,(q)=0, /gK}. (7.7) 

Next, we simply modify the two constrained updates in the above pseudo-code, given in Algorithm 1, to include 
the equality constraints. We then obtain a new generalized equality-inequality position update step (to replace 
line 5 in the pseudocode): 

q^{x: DiL^(q,x) + p + Ns(q)A + F(q)v=0, 0<A±gs(x)>0, f(x)=o} (7.8) 

and a corresponding momentum update step (to replace line 6 in the pseudocode): 

P^{y: y = Z)2i-./(q"",q) + Ns(q)Al + F(q)^, Ns(q)^M-V = 0, F(q)^M-V = o}. (7.9) 

Otherwise, the remaining pseudocode and all other previous discussion, derivation and guarantees (smooth- 
interval symplecticity, momentum conservation, etc.) continue to hold. 

8. Further Numerical Examples. To understand and examine the behavior of the algorithms proposed 
here, we investigate implementations of GVI over a range of additional numerical examples. Although a wide 
variety of predictors can be employed, in the following examples we exclusively apply the simple, uncon- 
strained, forward predictor, 

q'' = {x: DiL^(q',x) + p'=0}, (8.1) 



Here ■+ indicates the standard matrix pseudoinverse. 
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Fig. 8.1. Nonlinear Oscillator Example: In these figures we examine the behavior of various algorithms on the nonlinear oscillator 
example in two dimensions. All methods shown here begin with implicit midpoint as the base, unconstrained numerical method. See Figure 
\6. 1 \ (c} and §g. 0. J ^ or further details. Green lines and red lines plot energy and angular momentum respectively. Note that both energy and 
angular momentum should be preserved in this example. In (a) we plot the endpoint-constrained direct method. In (b) we plot the standard 
hybrid method. In (c) we plot the GVI integrator, where energy and angular momentum are both preserved long-term. Finally, in (A), to 
highlight the importance of enforcing nonlinear constraints appropriately, we plot the GVI integrator using linearized constraints. Note 
that, for this last example, the GVI method continues to preserve energy; however, linearized constraints will no longer preserve rotational 
symmetry and thus, correspondingly, angular momentum is no longer preser\>ed. 



and reserve an investigation of alternate predictors for later research. We implement implicit-midpoint and 
Stormer-Verlet based GVI-schemes using the implicit midpoint quadrature, 



L,,(q^q''+^)=/2L( 



q<r+l_l_qir qA:+i_q/; 



(8.2) 



and the Verlet quadrature, 



LM,^'^') = ^[L{o, 



, q'^+'-q^ 



-i(q 






)\ 



(8.3) 



respectively. 

8.0.1. Nonlinear Oscillator. To examine a nonintregable system and momentum conservation properties, 
we next consider the example of a nonlinear oscillator in M^. We start with two spheres in the plane, with 
respective positions qi and q2 in M^ and radii r\ and r2- The full configuration is then q = (qi,q2)^ with a 
corresponding nonlinear potential for the oscillator 



2 1^2 



y(q) = ^|lq,|l2(llq,ll'-l 

!=1 



(8.4) 



See the inset figure below. We then impose the hard-sphere, non-overlap constraint between the two spheres 
using the nonlinear, inequality constraint 



g(q)=!l qi-q2 II -n-'-2>o. 



(8.5) 



^(0, 



PiloT-^ 



Note that both the potential and constraint in this example are invariant with respect to rotations and thus 
angular momentum should be conserved. 

For this example all methods considered begin with implicit-midpoint as the base, 
unconstrained numerical method. We then set the step-size \.oh= 10^^ mass and sphere 
radii are set to m — 1 and r ~\, respectively, for both particles. The starting configura- 
tion is given such that the particles are initially separated at the non-equilibrium position, 
qj(0) = (0,-r-4x 10"')^, q2(0) = (0,r + 4x 10"')-^, while the initial momentum ini- 
tiates a counter-clockwise rotation with pi(0) == (1,0)"^ and P2(0) = (— 1,0)^. The system 
should then, as illustrated to the right, enter into a relative equilibrium composed of both 
periodic bouncing between the particles combined with a global rotation of the full system 
about the origin. 
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We compare the results of the endpoint-constrained direct method rl (i.e., imposing g(q'+^) > 0), the 
standard hybrid method, GVI appHed using the full nonlinear constraint given in (8.5 i, and GVI applied using 
a linearization of the constraints. We include this last example to highlight the potential pitfalls associated with 
applying constraint linearization since this a popular simplifying strategy in the contact mechanics literature. 
In Figure [8T| we plot the obtained energy in green and angular momentum in red for each method. 

In Figure |8Tf a), as might be expected, the endpoint-constrained direct method quickly dissipates all nor- 
mal modes of oscillation between the two spheres. After this initial dissipation, the spheres enter a relative 
equilibrium in which they orbit, in constant contact, with a purely smooth motion. This latter mode corre- 
sponds exactly to an endpoint-constrained direct method applied to a smooth system. As discussed in Section 



3.1 direct methods applied to smooth motion are likewise not symplectic and thus are still not guaranteed to 
preserve energy. This is reflected in the plot, where the system continues to generate correspondingly poor 
energy behavior. In particular, for this latter phase, we see slow but continual energy growth. 



In Figure 8.1 b) the standard collision integration method likewise displays poor energy behavior. As 



discussed in Section 1.3.2 energy drift (in the form of dissipation in this case) is generated, in collision inte- 
grators, by the interleaving of reflections with variable step-size, unconstrained integration steps. The reflection 
impulses, however, are applied along constraint gradients evaluated at the time of constraint activation and thus, 
despite the poor energetic behavior, momentum is conserved due to the rotational invariance of the constraint 
function. 



In Figures 8.1 c) and (d) we plot GVI using the true nonlinear constraints and linearized variants respec- 



tively. Our GVI method, plotted in Figure 8.1 c), enforces the true nonlinear constraints and correspondingly 



we note that both approximate energy and exact angular momentum are conserved long-term, while the com- 
bined nonsmooth, rotational and oscillatory modes are all preserved. Finally, we note that, when applied to 



linearized constraints, as in Figures 8.1 d), GVI continues to preserve energy but, due to constraint lineariza- 
tion, rotational invariance no longer holds and thus angular momentum is no longer conserved. This highlights 
the cost of linearizing constraints when applying geometric methods. 
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Fig. 8.2. We consider the combined smooth and nonsmooth trajectory of a constrained system composed of the same simple 
constrained mass-spring system considered in ^5.2\ and Figure [577] Unlike the previous example, however, we now initialize the first 
particle to an off sphere boundary position and set the second particle to an on-boundary position as illustrated in (a). As in ^5.2\ both 
particles are initialized with corresponding unit length momenta that point counter-clockwise and tangent to the boundary. This initiates 
a counter-clockwise rotation for both particle with an oscillatory motion between them via forces applied by the spring. The first particle, 
however, should now begin a nonsmooth oscillatory bouncing off the surface, while the trajectory of the second particle should remain 
smooth and on boundary throughout the simulation. Otherwise all other physical parameters remain unchanged from the last example. 
In (b) we plot the energy (in green) and the angular momentum (in red) of the direct-substitution, end-point constrained, method for the 
first 100 units of simulation time. We note the dissipation in both angular momentum and energy. This includes a rapid decay of all 
nonsmooth bouncing modes. In (c) we similarly plot the (approximate) preservation of energy (in green) and (exact) conservation of 
angular momentum (in red) obtained by GVI out to 1000 units of simulation time. See ^8. l\f'or further discussion. 



8.1. Smooth and Nonsmooth, Spring and Sphere. 



smooth and nonsmooth boundary modes, we reconsider our spring and sphere example from { 5.2 



As a first consideration of an example with combined 

In this 



example we leave all physical and time-step parameters unchanged from from S 5.2 except that we now initialize 
the example system to a combined smooth and nonsmooth, rotational trajectory. We initialize the first particle 
to an off sphere boundary position at qi(0) = (4,-1). We then set the second particle to an on boundary 
position, q2(0) = (3,-4). See Figure 8.2 a). As in E5.2 both particles are initialized with corresponding unit 



length momenta that point counter-clockwise and tangent to the boundary. This initiates a counter-clockwise 



*The midpoint-constrained direct method (i.e., imposing g( - — ^r^ ) ^ 0) was higlily unstable for this example. 



rotation for both particles with an oscillatory motion between them via forces applied by the spring. The first 
particle, however, should now begin a nonsmooth oscillatory bouncing off the surface, while the trajectory of 
the second particle should remain smooth and on boundary throughout the simulation. 



In Figure 8.2 (b) we plot the energy (in green) and the angular momentum (in red) of the direct, end-point 
constrained method for the first 100 units of simulation time. We note that decay in both angular momentum 
and energy for the direct method corresponds to the destruction of both the nonsmooth and smooth modes of 



the particles, as well as the rotational and oscillatory modes of the system. In Figure 8.2 (c) we similarly plot 



the (approximate) preservation of energy, in green, and the (exact) conservation of angular momentum, in red, 
obtained by our GVI method out to 1000 units of simulation time. Throughout this simulation, GVI maintains 
the nonsmooth bouncing mode of the first particle, the smooth-on-boundary mode of the second particle, as 
well as the rotational and oscillatory modes of the full system. 
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Fig. 8.3. Newton's Cradle Example: In (a) we compare the results of our impUcit-midpoint-based GVI method applied to a 
Newton's Cradle System over a range of decreasing step sizes, h = 3 x 10^^ (green), 2 x 10^^ (blue), 10^~ (red), and 5 x 10^' (black). 
We note that in the passage from h = 3x 10"^ to h= 10"^, we obtain the classic transition to stable energy behavior that is characteristic 
of symplectic methods. In (b) we illustrate the behavior of the integrator with a series of snapshots from the simulation. 

8.1.1. Newton's Cradle. Next we examine an example with combined smooth and nonsmooth boundary 
modes and the simultaneous enforcement of both equality and inequality constraints. We consider the classic, 
pendulum-based, Newton's Cradle system composed of n initially touching spheres subject to gravity. Each 
sphere / e [1,«] is constrained by a corresponding equality constraint 



f,(q)-||q,-p,||-Z=0, 



(8.6) 



to remain at constant length I from its hanging attachment point p,. Neighboring pairs /,y of spheres are then 
further constrained by non-overlap constraints 



g,j(q)=l|q,--q>ll -'•; - o > o. 



(8.7) 



Thus, for an n-sphere cradle we enforce n bilateral and n — 1 unilateral constraints and obtain a combination 
of smooth motion along the constraint boundary interfaces between neighboring spheres in moving and resting 
clusters, holonomic motion along the manifolds given by the pendulum constraints, and nonsmooth motion at 
each sphere impact. 

For this example we begin with implicit-midpoint as the base, unconstrained numerical method. We then 
set the masses, sphere radii, and pendulum lengths to m = 1, r = 0.25 and 1 — 1, respectively, for all spheres. 
We examine a five-ball cradle. The starting configuration is given such that all but the two leftmost spheres 
are at rest and in contact. The two leftmost spheres are initially pulled back into a non-equilibrium, contacting 
position, at an angle of —n/5 from rest. The system should then enter into the characteristic Newton's cradle 



behavior as illustrated by the simulation snapshots shown in Figure 8.3 (b) 



We compare the results of GVI over a range of decreasing step sizes. In Figure 8.3 (a) we plot the obtained 
energy for h — 3 x 10"^, 2 x 10"^, 10"^, and 5 x 10"^. We note that in the passage from h = 3 x 10"^ to 
h — 10^^ we obtain the classic transition to stable energy behavior that is characteristic of symplectic methods. 
Furthermore, as step size decreases, the cradle preserves the correct cyclical behavior, expected from Newton's 
cradle, for increasingly longer periods. Finally, we note that both the direct methods, and the hybrid schemes, 
along with generating their expected energy errors, do not predict the correct physical behavior of Newton's 
Cradle. 
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Fig. 8.4. Lennard-Jones Example: In these figures we examine the behavior of Stormer-Verlet-based GVI on a system of six 
particles in "S? connected sequentially by rods and interacting in a Lennard-Jones potential. In (a) we take a snapshot of the system 
configuration, (b) traces out the entire trajectory of all particle centers over 2000 units of simulation time, stepped ath= 10 , while (c) 
plots the change in energy over the same trajectory. 



8.1.2. Interacting Constrained Particles under Lennard-Jones. Finally, motivated by potential appli- 



cations in Molecular Dynamics |Frenkel and Smit|2001|, we consider the example of six spheres in M?, i.e.. 



(q{ , ..., qg )^ e M}^, connected sequentially by bilateral rod constraints, 

f/(q)==ll q,-q,+i 

Sphere-pairs are constrained by non-overlap conditions. 



/=0. 



g,/q)=||q,.-q,.||-2r>0, 



interact in a Lennard-Jones potential. 



V(q)-L4e 



12 



q,-q>lK Miq,--q,/ 

and are further constrained to lie in the interior of a bounding sphere of radius rg so that additionally. 



:>i.b 



(q) = -||q,.||-r + r,>0. 



(8.8) 



(8.9) 



(8.10) 



(8.11) 



In our simulation we employ reduced time units, (mc^/e)''^, with the mass of aU particles set to unity, o /2 
0.3, Vb 



r = U.:5, ri, — 5, and / = 1. Figure 8. LI (a) shows a snapshot of the system configuration. Figure 8. LI (b) 
traces out the entire trajectory computed by Storm er-Ver let-based GVI of all particle centers over 2000 units of 

^-2, 



simulation time, stepped at /i = 10 , while Figure 



8.1.1 



(c) plots the change in energy over the same trajectory. 



9. Discussion and Conclusion. We have introduced a fully nonsmooth, discrete Hamilton's Principle 
and an associated DELI integration formulation that enables the generation of Vl-based, geometric methods 
for the numerical integration of generalized inequality/equality-constrained, nonsmooth Hamiltonian systems. 
Adding additional structure to this framework we have obtained the first family of DELI-based integration 
methods: GVIs. By construction the GVI methods proposed here preserve momentum and equality constraints. 
They enforce simultaneous inequality constraints, obtain smooth unilateral motion along constraint boundaries 
and allow for both nonsmooth and smooth boundary approach and exit trajectories. Finally, we have shown 
that the proposed GVIs are symplectic over smooth-trajectory intervals and are observed to conserve energy 
approximately throughout. 

We have validated these properties in our numerical experiments and compared them against correspond- 
ing direct-substitution and collision-integration schemes in the literature. We have shown that our proposed 
methods resolve systems subject to multiple persistent, frequent and/or simultaneously active constraints. We 
have further tested these methods on difficult scenarios, where both smooth and nonsmooth active constraint 
modes occur simultaneously with high frequency; generally, in many of the examples presented above, during 
the majority of time steps. 

In many of the above examples, potentially difficult, nonlinear, numerical optimization problems were 
solved. In particular, we note that the complementarity conditions in our Discrete Smooth Integrator are, 
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unlike more standard KKT-type complementarity conditions, asymmetric with respect to constraint and con- 



straint gradient evaluations. As such, equations (|5.6|l - (5.9 1, do not correspond to the optimality conditions of 



standard nonlinear optimization problems | Bertsekas 1999| and instead require customized numerical solution 



strategies. In this work we have focused on investigating the challenges of inequality-constrained, numerical in- 
tegration and then on developing and validating our proposed methods. In future work we will explore efficient 
computational techniques for the scalable solution of these DELI-derived optimization problems. In practice, 
however, for the numerical experiments discussed here, we have found that relatively simple customized ap- 
proaches were possible that allowed for efficient solutions, while preliminary experiments investigating the 
possibility of extending these approaches to larger-scale problems are promising. 

Further, we look forward to examining the potential trade-offs between energy behavior and various de- 
grees of constraint enforcement guarantees. We also expect that a rigorous constraint error analysis, based 
on choice of predictor, should be a useful complement to the GVI methods proposed here. Finally, we note 
that GVI methods are the first of many possible instantiations of DELI methods. The geometric integration of 
inequality-constrained systems clearly remains a challenging area of investigation and we hope that the DELI 
framework will lead to further fruitful explorations, and the development of new complementary numerical 
methods. 
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Appendix A. The Discrete Generalized Reflection Operator. 



For an arbitrary choice of energy function, ^(q, p), we invoke the variational model from Kaufman et al. 



|2010|, in which a generalized reflection operator is interpreted as a sequence of energy preserving, constraint 



satisfying projections. Following their time-continuous case derivation we then generate the corresponding 
discrete generalized reflection analogue as the discrete, energy preserving, constraint satisfying projection se- 
quence defined in Algorithm 2 below. 

Algorithm 2 Discrete_Generalized_Ref lection(q,p,K,£') 



while true do 

for ^ in K do 

ifVg^.(q)^M"'p<0 then 

end if 
end for 

if V 7^ ©then 

A^{y: £(q,p + Gv(q)y)=£(q,p), Gv(q)^M-' (p + Gv(q)y) > 0, y > o} 

p^ p + Gv(q)A 
else 

return p 
end if 
end while 



Theorem K.\. At termination. Algorithm 2 satisfies all discrete jump conditions. 



Proof By construction, each iteration satisfies conditions (6.8 i and (6.9 1. Termination guarantees condi- 
tion ( |6J] l . D 

When the energy function is given by the separable continuous Hamiltonian, e.g., E = H, the discrete 
generalized reflection operator is identical to the time-continuous case considered in Kaufman et al. |2010J. In 



this case each solve of line 9, in the above pseudo-code, is obtained by the non-negative, least-squares solution 
of the normal equation: 

Gv(q)^M-'Gv(q)A = -2Gv(q)^M-ip, A > 0. (A.l) 

More generally, whenever the Energy function is quadratic in p, e.g., as in the Stormer-Verlet's second-order 
numerical Hamiltonian | Bond and Leimkuhler 2007], the solution to line 9 is likewise given by a non-negative. 



linear least-squares system. Finally, for possible energy functions containing higher order polynomials in p, 
more complex, variational treatments are required. 
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Finally, we note that, as discussed above, GVI is agnostic to the choice of reflection model as long as 



a well-posed solution exists that satisfies the discrete jump conditions, (6.7 i - (6.9 1. As in the case of the 
Generalized Reflection operator, alternate existing multi-impact models can be similarly discretized to do so. 
As a particularly compact example, Moreau's elastic multi-impact solution [ Moreau|1988| can be discretized 
to obtain the energy preserving discrete reflection given in Algorithm 3 below. 



Algorithm 3 Discrete_Moreau_Ref lection(q,p,K,£') 



1: A^{y: £(q,p + GK(q)y)=£(q,p), GK(q)^M-i (p + GK(q)y) > 0, y > o} 
2: p^p + GK(q)A 
3: return p 



Appendix B. Direct Substitution Schemes. 

As discussed above, an extended value potential for inequality constrained systems can be composed by 
concatenating the unconstrained Hamiltonian system's potential with the extended value indicator function. 
Following this observation, the direct-substitution approach in symplectic methods [ Kane et al.|1999 Stewart 
2000 Pandolfi et al. 12002) [Deuflhard et al.|[2008j | is to compose direct extensions by treating the nonsmooth 
constraint force, — (9/^(q), as simply an additional force to be handled by each numerical method in the usual 
manner. 

B.l. Direct Integration Examples. To illustrate this approach, in the following, we will present examples 
using two popular integration approaches: the implicit midpoint method and the Newmark family of integration 
algorithms. 

B.1.1. Implicit Midpoint Method. The implicit midpoint method is given by the discrete, momentum, 
phase-space pair 



q'+l=q' + jM-l(p'+'+p') 



p'+i = p'_/,vy( 



q'+'+q' 



(B.l) 
(B.2) 



A direct-substitution, nonsmooth-constrained extension to the midpoint rule is generated by adding the 
extended value indicator to the unconstrained system's potential. Applying a generalized gradient, this gives 



q'+i=q' + _M-'(p'+'+pO, 

„t+\ , t r+l I r 

p'+'ep'-wy(^^)-/.5/A(^^^). 



(B.3) 
(B.4) 



An alternate, end-point constrained variant of implicit midpoint can be obtained by simply evaluating 
constraint forces at the end of the time-step [Stewart 2000J . This leads to the modified midpoint method 
momentum update 



p'+i e p'-wy(5^^^)-/za/A(q'+'] 



(B.5) 



B.1.2. Newmark Methods. The Newmark family of integration schemes are given by the discrete, ve- 
locity, phase-space pair 



/z2 

2 



q'+i=q'+/iq'-— (l-2j3) M-^VV[c()+2^ U-^VV{a!+^) 



q'+^=q'-/2 (1-7) M-iVy(q') + 7M-iVy(q'+i) 
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(B.6) 
(B.7) 



Recall that Newmark methods vary j3 G [0, 5] and 7 e [0, 1]. Newmark is second order accurate if and only if 
7=5, otherwise it is only consistent. Implicit Newmark with jS = | and Y—j gives the implicit trapazoidal 
method, while fi —0 and 7=5 gives explicit Newmark. 

A direct-substitution, nonsmooth-constrained extension of the Newmark Family is generated by adding the 
extended value indicator to the unconstrained system's potential. Applying a generalized gradient, this gives 



q'+'eq'+/iq'--[(l-2j3) M"' (vy(q') +<9/A(q')) +2)3 M'' (vviq'+^) +dl^{q'+^) 
q'+'eq'-h\{l-r) M-i(vy(q') + <9/A(q'))+7M-'(vy(q'+i) + <9/A(q'+') 



(B.8) 
(B.9) 



An Implicit/Explicit (IMEX) variant of this extension was developed by Kane et al. 1 1999 1 in which both 



/3 and 7 are enforced as fully implicit values (i.e., j3 = 5 , 7 = 1) for just the nonsmooth portion of the composite 
potential. This generates a family of nonsmooth IMEX Newmark integrators given by 



U2 tp. 

q'+^eq'+/!q' (l -2j3) M"' Vy(q') +2j3 M"' VV (q'+i) U-'^dI/^{a!+^), (B.IO) 



q'+^eq' 



-h 



(1-7) M-'Vy(q') + 7M-iVy(q'+i) -hU-^dIi^{q'+'') 



(B.ll) 



B.2. Variational Interpretation of Direct Nonsmooth Integrators: Optimization Forms. Given the 
explicitly variational definition of the generalized gradient [Rockafellar and Wets 1998 J , direct, nonsmooth 
constrained extensions of symplectic methods, as described above, generally reduce to constrained, nonlinear 
minimizations of the form 



q'+l=argmin{e(q',q,/i): f (q',q,/z) € A J. 



(B.12) 



Interpreted in this equivalent form, symplectic methods, applied to Hamiltonian systems, generate an 
objective function, e, generally polynomial in q, and a vector-valued, constrained configuration function, f, 
linear in q. 

B.2.1. Implicit Midpoint Example. As an example, consider the nonsmooth midpoint rule. We substi- 



tute Equation ( |B.4| i into \B3\ to obtain the DI 

/,2 

q'+' eq'+/iM-'p'-— M-'Vy( 



.,?+! 



V+1 



-)• 



(B.13) 



2 ' 2 ' 2 '^^ 2 

Applying the definition of the generalized gradient, we then obtain an equivalent local, nonlinear minimization 



that corresponds to the template given by (B.12 1, 
q'^ = argmin-^ -q^Mq — q^(q' 



-h^vC- 



-hM-^p') 



Rearranging ( |B.3[ ) the corresponding midpoint method, momentum update is then given by 

2. 



^):(5^)eA}. (B.14) 



^r+l 



M(q'+l-q')-p'. 



(B.15) 



Newmark and other direct methods follow similarly. 



B.3. Direct Method Behavior. Unfortunately, many of the characteristic advantages of symplectic meth- 
ods are lost in the above direct nonsmooth treatment. The chief observation is, in the smooth setting, the 
good behavior of symplectic methods is ascribed, via backwards error analysis, to the existence of trajectory 
shadowing numerical Hamiltonians JHairer et al.|2002] . More specifically, the flow of fixed-step, symplectic 
methods follow, up to truncation, a numerical Hamiltonian system that is constructed via a series expansion. 
Because the numerical Hamiltonian is obtained under a smoothness assumption, which does not hold in the 
nonsmooth case, standard backwards error analysis guarantees no longer hold for the inequality constrained 
systems | |Stewart|2000l|Bond and Leimkuhler|2007 1. 
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Despite these issues, during free motion (i.e., when constraints are inactive), these nonsmooth constrained 
symplectic methods trivially reduce to their corresponding unconstrained method and are thus still guaranteed 
to shadow numerical Hamiltonians. At each constraint boundary event, however, this guarantee is lost, and the 
solution obtained no longer matches the shadow trajectory. 

In practice, it is exactly at these boundary events that such direct methods obtain a wide variety of incorrect 
and undesirable behaviors. Constraint enforcement generates dissipative trajectories for some direct methods, 
energy growth in others, and in many cases produces effectively nondeterministic restitution behavior both 
when varying time-step sizes across simulations and within a single fixed step-size simulation |Stewart|2000|. 



B.4. No Energy Conservation. Considering (B.12i we, in particular, note that, perhaps surprisingly, a 



large reason for these difficulties stems directly from the underlying variational structure of these methods. 

Theorem B.l. All direct, one-step, nonsmooth symplectic methods that can be posed in the form of 
\B.12\ , can not, in general, preserve, either approximately, nor exactly, the energy of the underlying uncon- 
strained symplectic method. 

Proof. For all such methods, ( |B.12| l can be rewritten, with respect to the constraint functions, as 

q'+'=argmin{e(q',q,/2): g((l - a)q' + aq) >o}, (B.16) 

for some fixed a £ (0, 1]. 

We then note that we retrieve the solution of the underlying, unconstrained symplectic method by solving 
for q in 

V^e{q',q,h)=0. (B.17) 

The constrained variant will return this solution for any step in which the unconstrained solution returns a non- 
violating c onfiguration. M ore generally, formulating the Karush-Kuhn-Tucker (KKT) first order optimality 



conditions |Bertsekas 1999J, for the above minimization, we obtain, for optimal q'+' = q, 

V^e{q',q,h) - aVg((l - a)q' + aq^X = 0, (B.18) 

0<A_Lgf(l-a)q' + aq) >0. (B.19) 



We then consider any time step in which an unconstrained solution violates constraints. The above KKT 
system implies that constraint forces are appUed along constraint gradients and, more specifically, are of the 
form a V g A , A > 0. 

The complementarity term of the above KKT system then additionally requires, component-wise, that a 
constraint force magnitude, applied along the gradient Vgi, can be nonzero (and thus enforce the corresponding 
constraint, gf) if and only if 

g,((l-a)q' + aq'+')-0. (B.20) 

Thus constraining forces can only be applied along constraint gradients (to enforce feasibility) if the position 
at which the constraint force is evaluated, (1 — a)q' + aq'+', lies exactly on the corresponding constraint's 
boundary. Since these methods demand constraint enforcement by construction, they implicitly enforce this 
additional condition as a secondary constraint. In turn, this implicit constraint overrides the underlying sym- 
plectic structure of the base, unconstrained method. In particular, the magnitude of the corresponding constraint 
force, that achieves this secondary constraint, scales independently of time-step size and energy. Instead, these 
constraint force magnitudes scale with constraint geometry and state, so that an arbitrarily large perturbation 
may be applied to the discrete system's energy and momentum in order to place the configuration at which 
constraint forces are evaluated, (1 — a)q' + aq'+^ on the admissible set boundary. D 

Appendix C. The Discrete-Smooth Integrator. We now show symplecticity and momentum conservation 
for the Discrete-Smooth Integrator (DSI) by proving Theorem 5.1. 

Proof. As a preliminary, starting at time f+, we first define an Oracle Set composed of constraint indices 

that will be active at the end of the current time step, 0(q'+') = {/ : giiq''^^) = 0}. The intersection of the 

oraclesetwiththe5moof/!5ef,§(q',p'^), defined in (6.1), gives I ^'^ 0(q'+') n §(q',p'^). 
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Now consider the optimal A and ji that obtain (q'+^, p'+') in (6.3) - (6.6). We then observe that setting 
7= A^ and 5 = /i we also satisfy the following equality constrained integration system 

DiL,,(qSq'+i) + p' + Gi(q')7= 0, (C.l) 

gl(q'+')-0, (C.2) 

p'+l=D24/(q',q'+') + Gi(q'+')5, (C.3) 

Gi(q'+i)^M-ip'+'=0. (C.4) 

In particular, plugging 7 == A ^ and 8 = jX into the above system gives exactly the same (q'+ ' , p'+ ' ) solution as 
the DSI system. 

Next we observe that Equations (|5.6[) through (|5.9[) give the position-momentum form of the equality con- 



f |M^ 



strained discrete Hamiltonian map of Marsden and West 1 2001 (3.5.2a) - (3. 5. 2d)]. Symplecticity and momen- 



tum conservation then follow directly by noting that restricting DSI steps to constraints indexed in §(q', p' ) 
guarantees initial position and momentum feasibility ]Marsden and West|2001| (3.5.1)] for the corresponding 
equality constrained discrete Hamiltonian maps, at all steps. D 
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